#---
#FUNCTIONS: ESTIMATE INTENSITY FUNCTION USING LOCAL LINEAR LOCAL-EM ALGORITHM
#---

#quad5pt: computes Gaussian quadrature for 5 points
quad5pt <- function(J)
{
	JintervalWidth <- diff(c(0, J))
	nJ <- length(J)

	absc <- c(-0.906179845938664,
        -0.538469310105683,
		0.000000000000000,
        0.538469310105683,
        0.906179845938664)
	wgt <- c(0.236926885056189,
        0.478628670499366,
        0.568888888888889,
        0.478628670499366,
        0.236926885056189)

	nQuadpts <- length(absc)
	halfw <- 0.5 * JintervalWidth
	midpts <- 0.5 * (J + c(0, J[-nJ]))
	x <- rep(midpts, each = nQuadpts) + rep(halfw, each = nQuadpts) * rep(absc, nJ)
	return(list(nQuadpts = nQuadpts, lattice = x, wgt = wgt))
}

#quad10pt: computes Gaussian quadrature for 10 points
quad10pt <- function(J)
{
	JintervalWidth <- diff(c(0, J))
	nJ <- length(J)

	absc <- c(-0.973906528517172,
        -0.865063366688985,
        -0.679409568299024,
        -0.433395394129247,
        -0.148874338981631,
        0.148874338981631,
        0.433395394129247,
        0.679409568299024,
        0.865063366688985,
        0.973906528517172)
	wgt <- c(0.066671344308688,
        0.149451349150581,
        0.219086362515982,
        0.269266719309996,
        0.295524224714753,
        0.295524224714753,
        0.269266719309996,
        0.219086362515982,
        0.149451349150581,
        0.066671344308688)

	nQuadpts <- length(absc)
	halfw <- 0.5 * JintervalWidth
	midpts <- 0.5 * (J + c(0, J[-nJ]))
	x <- rep(midpts, each = nQuadpts) + rep(halfw, each = nQuadpts) * rep(absc, nJ)
	return(list(nQuadpts = nQuadpts, lattice = x, wgt = wgt))
}

#moment generating function for Gaussian distribution
iMGF.gaussian <- function(ub, t = 0)
{
	exp(0.5 * t^2) * pnorm(ub, mean=t)
}

iMGF.gaussian.first.derivative <- function(ub, t = 0)
{
	exp(0.5 * t^2) * (t * pnorm(ub, mean = t) - dnorm(ub, mean = t))
}

iMGF.gaussian.second.derivative <- function(ub, t = 0)
{
	exp(0.5 * t^2) * ((1 + t^2) * pnorm(ub, mean = t) - (ub + t) * dnorm(ub, mean = t))
}

#score function
score.fun <- function(a0, a1, lattice, bandwidth, L, R, em.estimate, entry1, entry2)
{
	if(any(L > R))
	{
		stop("Right endpoints have to be greater than left ones.")
	}

	if(length(R) == length(L))
	{
		nrows <- length(R)
		hessian.mat <- diag(0, 2 * length(lattice))
		indexA0 <- 2 * (1:length(lattice)) - 1
		indexA1 <- 2 * (1:length(lattice))
	}

	else
	{
		stop("Lengths of left- and right-end points are different.")
	}

	if(length(a0) != length(lattice))
	{
		stop("Constant coefficient and grid points are of different length.")
	}
	
	if(length(a1) != length(lattice))
	{
		stop("Linear coefficient and grid points are of different length.")
	}
	
	#first entry of the score function
	firstTerm <- outer(R, lattice, function(x, y) iMGF.gaussian(ub = (x - y)/bandwidth)) -
		outer(L, lattice, function(x, y) iMGF.gaussian(ub = (x - y)/bandwidth))
	firstTerm <- as.vector(as.vector(em.estimate) %*% diag(1/(R - L)) %*% firstTerm) + as.vector(entry1)
	
	#second entry of the score function
	secondTerm <- outer(R, lattice, function(x, y) iMGF.gaussian.first.derivative(ub = (x - y)/bandwidth)) -
		outer(L, lattice, function(x, y) iMGF.gaussian.first.derivative(ub = (x - y)/bandwidth))
	secondTerm <- as.vector(as.vector(em.estimate) %*% diag(bandwidth/(R - L)) %*% secondTerm) + as.vector(entry2)
	
	return(as.vector(rbind(firstTerm, secondTerm)))
}

#hessian function
hessian.fun <- function(a0, a1, lattice, bandwidth, L, R, offsets)
{
	if(length(R) == length(L))
	{
		nrows <- length(R)
		hessian.mat <- diag(0, 2 * length(lattice))
		indexA0 <- 2 * (1:length(lattice)) - 1
		indexA1 <- 2 * (1:length(lattice))
	}

	else
	{
		stop("Lengths of left- and right-end points are different.")
	}
	
	if(length(a0) != length(lattice))
	{
		stop("Constant coefficient and grid points are of different length.")
	}

	if(length(a1) != length(lattice))
	{
		stop("Linear coefficient and grid points are of different length.")
	}

	tmpA1 <- rep(a1, each = nrows) * bandwidth
	
	#second derivative w.r.t a0
	tmp11 <- outer(R, lattice, function(x, y) iMGF.gaussian(ub = (x - y)/bandwidth, t = tmpA1)) -
		outer(L, lattice, function(x, y) iMGF.gaussian(ub = (x - y)/bandwidth, t = tmpA1))
	tmp11 <- -exp(a0) * as.vector(offsets %*% tmp11)
	diag(hessian.mat)[indexA0] <- tmp11
	
	#cross derivatives
	tmp12 <- outer(R, lattice, function(x, y) iMGF.gaussian.first.derivative(ub = (x - y)/bandwidth, t = tmpA1)) -
		outer(L, lattice, function(x, y) iMGF.gaussian.first.derivative(ub = (x - y)/bandwidth, t = tmpA1))
	tmp12 <- -bandwidth * exp(a0) * as.vector(offsets %*% tmp12)
	hessian.mat[(row(hessian.mat) == col(hessian.mat) + 1) & (row(hessian.mat) %% 2 == 0)]<- tmp12
	hessian.mat[(row(hessian.mat) == col(hessian.mat) - 1) & (row(hessian.mat) %% 2 == 1)]<- tmp12
		
	#second derivative w.r.t a1
	tmp22 <- outer(R, lattice, function(x, y) iMGF.gaussian.second.derivative(ub = (x - y)/bandwidth, t = tmpA1)) -
		outer(L, lattice, function(x, y) iMGF.gaussian.second.derivative(ub = (x - y)/bandwidth, t = tmpA1))
	tmp22 <- -(bandwidth^2) * exp(a0) * as.vector(offsets %*% tmp22)
	diag(hessian.mat)[indexA1] <- tmp22
	
	if (det(hessian.mat) != Inf)
	{
		hessian.inv <- solve(hessian.mat)
	}

	else
	{
		hessian.inv <- matrix(0, nrow=nrow(hessian.mat), ncol=ncol(hessian.mat))
		for(k in 0:(length(R) - 1))
		{
			tmpInd <- (20*k+1):(20*k+20)
			hessian.inv[tmpInd, tmpInd]<- solve(hessian.mat[tmpInd, tmpInd])
		}   
	}

	#tmp11 and tmp 12 are needed to calculate the score function
	return(list(hessian = hessian.mat, hessian.inv = hessian.inv, entry1 = tmp11, entry2 = tmp12, entry3 = tmp22))
}

#ems.constructor.linear: computes smoothing matrix
#panel.data - panel data
#bandwidth - bandwidth
#quad.order - number of points
ems.constructor.linear <- function(J, lattice, bandwidth, a1, offsets, quadwgts)
{
	nJinterval <- length(J)
	JintervalWidth <- diff(c(0,J))
	tmpA1 <- rep(a1, each = nJinterval) * bandwidth
	
	#construct the smoothing matrix at different x (for the final intensity estimate)
	num <- outer(J, lattice, function(x, y) iMGF.gaussian(ub = (x - y)/bandwidth)) - 
		outer(c(0,J[-nJinterval]), lattice, function(x, y) iMGF.gaussian(ub = (x - y)/bandwidth))
	den <- outer(J, lattice, function(x, y) iMGF.gaussian(ub = (x - y)/bandwidth, t = tmpA1)) -
		outer(c(0,J[-nJinterval]), lattice, function(x, y) iMGF.gaussian(ub = (x - y)/bandwidth, t = tmpA1))
	K_x.mat <- diag(offsets/JintervalWidth) %*% num %*% diag(1/as.vector(offsets %*% den))
	
	#construct the smoothing matrix used in the proposed EMS algorithm by Gaussian quadrature rule
	K.mat <- K_x.mat %*% kronecker(diag(0.5 * JintervalWidth, nrow = nJinterval), quadwgts)

	return(list(K_x.mat = K_x.mat, K.mat = K.mat))
}

#local.linear.nhpp: iteratively compute intensity function estimates for nonhomogeneous Poisson process
#panel.data - panel data
#bw - bandwidth
#quad.order - number of points
local.linear.nhpp <- function(panel.data, bw, tol = 1.0e-8, max.iter = 100, quad.order = NULL, partition = NULL, Lambda.ini = NULL)
{
	if(is.null(partition))
	{
		J <- sort(union(panel.data$L, panel.data$R))
	}
	
	else
	{
		J <- partition
	}
	
	if(any(J == 0))
	{
		J <- setdiff(J, 0)
	}
	
	nJinterval <- length(J)
	JintervalWidth <- diff(c(0,J))

	if(quad.order == 10)
	{
		quadRule <- quad10pt(J)
	}
	
	else
	{
		quadRule <- quad5pt(J)
	}

	#construct the incident matrix
	censor.mat <- outer(tapply(panel.data$CT, panel.data$ID, unique), J, ">=")
	offsets <- apply(censor.mat, 2, sum)
	tmp.data<- panel.data[which(panel.data$M!=0),]
	incident.mat<- outer(tmp.data$R, J, ">=") - outer(tmp.data$L, J, ">=");
	
	if(is.null(Lambda.ini))
	{
		Lambda.old <- tmp.data$M %*% incident.mat %*% diag(1/offsets)
	}
	
	else
	{
		Lambda.old <- t(as.matrix(Lambda.ini))
	}
	
	#EM iteration
	a0.old <- kronecker(log(2 + as.vector(Lambda.old)), rep(1, quadRule$nQuadpts))
	a1.old <- rep(0, length(a0.old))
	a.old <- as.vector(rbind(a0.old, a1.old))
	a0 <- rep(0, length(a0.old))
	a1 <- rep(0, length(a0.old))
	a <- as.vector(rbind(a0, a1))
	Lambda <- matrix(0, nrow = 1, ncol = nJinterval)
	niter <- 1
	
	repeat 
	{
		V.inv <- diag(1/as.vector(Lambda.old %*% t(incident.mat)))
		tmpEM <- tmp.data$M %*% V.inv %*% incident.mat %*% diag(as.vector(Lambda.old))
	  
		#Newton Raphson algorithm
		nr.iter <- 0
		repeat
		{
			a0.old <- a.old[seq(1, length(a), by = 2)]
			a1.old <- a.old[seq(2, length(a), by = 2)]
			hessian <- hessian.fun(a0 = a0.old, a1 = a1.old, lattice = quadRule$lattice, bandwidth = bw, 
				L = c(0,J[-length(J)]), R = J, offsets = offsets)
			score <- score.fun(a0 = a0.old, a1 = a1.old, lattice = quadRule$lattice, bandwidth = bw, 
				L = c(0,J[-length(J)]), R = J, em.estimate = tmpEM, entry1 = hessian$entry1, entry2 = hessian$entry2)
			a <- a.old - as.vector(hessian$hessian.inv %*% score)

			if((sum((a.old - a)^2) > tol) & (nr.iter <= max.iter))
			{
				a.old <- a
				nr.iter <- nr.iter + 1
			}
	
			else
			{
				a0 <- a[seq(1, length(a), by = 2)]
				a1 <- a[seq(2, length(a), by = 2)]
				break
			}
		}
		
		#construct the smoothing matrix
		ems.input <- ems.constructor.linear(J = J, lattice = quadRule$lattice, bandwidth = bw, a1 = a1, 
			offsets = offsets, quadwgts = quadRule$wgt)
		
		#smoothing step
		Lambda <- (tmpEM %*% diag(1/offsets)) %*% ems.input$K.mat
		if(sum((Lambda.old - Lambda)^2) > tol)
		{
			Lambda.old <- Lambda
			niter <- 1 + niter
		}
	
		else
		{
			lambda <- (tmpEM %*% diag(1/offsets)) %*% ems.input$K_x.mat
			break
		}
	}
	
	return(list(Lambda.est = as.vector(Lambda), lambda.est = as.vector(lambda), nIter = niter, lattice = quadRule$lattice))
}
